!  FUNCTIONS/SUBROUTINES exported from SolDensMatrixFinalRound.dll:
module CXXWRAPPER
use iso_c_binding, only: c_double, c_int
implicit none
	interface
		!subroutine cxxgemv(matrix, x, y, m) bind(C, name="cxxgemv")
		!	import :: c_double, c_int
		!	real(c_double), intent(in), dimension(*) :: matrix
		!	real(c_double), intent(in), dimension(*) :: x
		!	real(c_double), intent(inout), dimension(*) :: y
		!	integer(c_int), value :: m
		!end subroutine cxxgemv

		!subroutine cxxgemv(y, matrix, x, n, lda) bind(C, name="cxxgemv")
		!	import :: c_double, c_int
		!	real(c_double), intent(inout), dimension(*) :: y
		!	real(c_double), intent(in), dimension(*) :: matrix
		!	real(c_double), intent(in), dimension(*) :: x
		!	integer(c_int), value :: n
		!	integer(c_int), value :: lda
		!end subroutine cxxgemv

		subroutine cxxgemv(y, matrix, x_full, x_start, n, lda) bind(C, name="cxxgemv")
			import :: c_double, c_int
			real(c_double), intent(inout), dimension(*) :: y
			real(c_double), intent(in), dimension(*) :: matrix
			real(c_double), intent(in), dimension(*) :: x_full
			integer(c_int), value :: x_start
			integer(c_int), value :: n
			integer(c_int), value :: lda
		end subroutine cxxgemv

		subroutine cxxgetrf(matrix, m) bind(C, name="cxxgetrf")
			import :: c_double, c_int
			real(c_double), intent(inout), dimension(*) :: matrix
			integer(c_int), value :: m
		end subroutine cxxgetrf

		subroutine cxxgetrs(matrix, b, m) bind(C, name="cxxgetrs")
			import :: c_double, c_int
			real(c_double), intent(in), dimension(*) :: matrix
			real(c_double), intent(inout),  dimension(*) :: b
			integer(c_int), value :: m
		end subroutine cxxgetrs

		!subroutine cxxgesv(matrix, ipiv, rhs, m) bind(C, name="cxxgesv")
		!	import :: c_double, c_int
		!	real(c_double), intent(inout), dimension(*) :: matrix
		!	integer(c_int), intent(out),  dimension(*) :: ipiv
		!	real(c_double), intent(inout), dimension(*) :: rhs
		!	integer(c_int), value :: m
		!end subroutine cxxgesv

	end interface
end module

	!subroutine mulMatrixVector(matrixSize,a,x,y)
	!	use CXXWRAPPER
    !!矩阵向量乘
	!	! use BLAS95
	!!   Expose subroutine mulMatrixVector to users of this DLL
	!  !DEC$ ATTRIBUTES DLLEXPORT::mulMatrixVector
	!	implicit none
	!	integer, intent(in) :: matrixSize
	!	real(kind=8), intent(in) :: a(matrixSize,matrixSize)
	!	real(kind=8),  intent(in) :: x(matrixSize)
	!	real(kind=8),  intent(inout) :: y(matrixSize)
	!
	!	! call gemv(a,x,y)
	!	call cxxgemv(a, x, y, matrixSize)
	!
	!end subroutine mulMatrixVector
	
	subroutine mulMatrixVector(matrixSize,a,matrixCapacity,x,x_start,x_end,y)
		use CXXWRAPPER
		!DEC$ ATTRIBUTES DLLEXPORT::mulMatrixVector
		implicit none
		integer, intent(in) :: matrixSize !矩阵尺寸
		integer, intent(in) :: matrixCapacity !储存矩阵的二维数组的尺寸
		real(kind=8), intent(in) :: a(matrixCapacity,matrixCapacity)!储存矩阵的二维数组
		real(kind=8), intent(in) :: x(:)!储存向量的一维数组
		integer, intent(in) :: x_start!向量在一维数组中所占空间的起始index
		integer, intent(in) :: x_end!向量在一维数组中所占空间的结束index
		real(kind=8), intent(inout) :: y(matrixSize)!输出的矩阵向量乘结果
		!call gemv(a(1:matrixSize,1:matrixSize),x(x_start:x_end),y)!上传作品需要替代的部分
		!call cxxgemv(y, a, x(x_start:x_end), matrixSize, matrixCapacity)
		call cxxgemv(y, a, x, x_start - 1, matrixSize, matrixCapacity) ! 1 indexing
		!call cxxgemv(y, a, x(x_start:x_end), 0, matrixSize, matrixCapacity) ! 1 indexing

	end subroutine mulMatrixVector

	SUBROUTINE getLU(matrixSize,a,ipiv)
		use CXXWRAPPER
    !LU分解
    
		! use LAPACK95
	
	  !DEC$ ATTRIBUTES DLLEXPORT::getLU
	
      IMPLICIT NONE
		integer, intent(in) :: matrixSize
        REAL(KIND=8),INTENT(INOUT)    :: a(matrixSize,matrixSize)
        INTEGER,intent(out)           :: ipiv(matrixSize)
		  
        !INTEGER                       :: INFO
          
		  
		  call cxxgetrf(a, matrixSize)
          !调用LAPACK函数进行矩阵LU分解, 将分解结果存在a里面
        !   CALL GETRF(a,IPIV,INFO)                      
          !IF(INFO/=0) THEN
          !    WRITE(14,'(A12,I5)') "DGETRF INFO=",INFO
          !    STOP
          !ENDIF
	END SUBROUTINE getLU
	
	
	SUBROUTINE forwardBackWardSubstitute(matrixSize,a,ipiv,b)
		use CXXWRAPPER
    !LU分解之后的向前向后替代
    
		! use LAPACK95
		!DEC$ ATTRIBUTES DLLEXPORT::forwardBackWardSubstitute
	
      IMPLICIT NONE
		integer, intent(in) :: matrixSize
        REAL(KIND=8),intent(in)    :: a(matrixSize,matrixSize)
        INTEGER,intent(in)          :: ipiv(matrixSize)
		REAL(KIND=8),intent(inout) ::  b(matrixSize)!输入右端项，输出求解结果
		  

		call cxxgetrs(a, b, matrixSize)
		! call getrs(a,ipiv,b)
	
    end SUBROUTINE forwardBackWardSubstitute
	
	

    SUBROUTINE solveLinearEquations(matrixSize,a,b)
		use CXXWRAPPER
    !求解线性方程组
    
		!use LAPACK95
		!DEC$ ATTRIBUTES DLLEXPORT::solveLinearEquations
	
      IMPLICIT NONE
		integer, intent(in) :: matrixSize
        REAL(KIND=8),INTENT(IN)    :: a(matrixSize,matrixSize)
		REAL(KIND=8),INTENT(INOUT) ::  b(matrixSize)!输入右端项，输出求解结果

		!REAL(KIND=8) :: aCopy(matrixSize,matrixSize)
		REAL(KIND=8),allocatable :: aCopy(:,:)
		integer :: ipiv(matrixSize)

		ALLOCATE(aCopy(matrixSize, matrixSize))

		aCopy = a

		print*, "gesv should not be called"
		!call cxxgesv(aCopy, ipiv, b, matrixSize)
		!call gesv(aCopy, b)

		DEALLOCATE(aCopy)
	
	end SUBROUTINE solveLinearEquations
	

